library(fmsb)
library(magick)
library(ggradar)
library(ggplot2)
library(cluster)
library(grid)
library(reshape2)  # for melt function
library(pheatmap)
library(RColorBrewer)

workDir <- c("/Users/meijian/Work/SIMPLE20181102/parameters_radarplot")
setwd(workDir)
speciesPara <- as.data.frame(t(read.csv("/Users/meijian/Work/SIMPLE20181102/VACS-SIMPLE_calibrated_parameters_for_all_crops.csv", 
                                        header=TRUE, sep=",", stringsAsFactors=FALSE, row.names = 1)))
paraDef <- c("Thermal time requirement from sowing to maturity", "Harvest index", "Thermal time requirement after sowing\nfraction of light interception to reach 50%", 
             "Thermal time requirement from maturity backwards\nfor light interception to reach 50%", "Base temperature", "Optimal temperature", 
             "Radiation use efficiency", "The maximum daily increase in I50B due to heat stress", "The maximum daily increase in I50B due to water stress",
             "Threshold for daily Tmax to start\naccelerating senescence due to heat stress", "Daily Tmax threshold when RUE becomes 0 due to heat stress",
             "Relative increase in RUE per 100ppm elevated CO2 above 350ppm", "Sensitivity of RUE to drought stress")
paraDef <- expression(T["sum"], HI, I["50A"], I["50B"], T["base"], T["opt"], RUE, I["50maxH"], I["50maxW"], T["max"], T["ext"], S["CO2"], S["water"]) # to create subscription
#typeNum <- c(6, 2, 7, 2, 5, 4)
typeNum <- c(6, 7, 2, 5, 4) # no fruits
donut_data <- data.frame(
  # category=rev(c("Cereals", "Fruits", "Legumes", "Oilseeds", "Roots/Tubers", "Vegetables")),
  category=rev(c("Cereals", "Legumes", "Oilseeds", "Roots/Tubers", "Vegetables")), #  no fruits
  count=rev(typeNum)
)
# Compute percentages
donut_data$fraction <- donut_data$count / sum(donut_data$count)
# Compute the cumulative percentages (top of each rectangle)
donut_data$ymax <- cumsum(donut_data$fraction)
# Compute the bottom of each rectangle
donut_data$ymin <- c(0, head(donut_data$ymax, n=-1))
# text location in y axis for crop categories
donut_data$label_position <- (donut_data$ymax + donut_data$ymin) / 2
# Set the angle for text labels, this is in reverse order of crop types!
# donut_data$mid_angle <- c(325,265,215,150,90,30)
donut_data$mid_angle <- c(325,255,200,135,35) # no fruits

## radar/spider plot
create_beautiful_radarchart <- function(data, color = "#00AFBB", 
                                        vlabels = colnames(data), vlcex = 0.5,
                                        caxislabels = NULL, title = NULL, ...){
  radarchart(
    data, axistype = 1,
    # Customize the polygon
    # pcol = color, pfcol = scales::alpha(color, c(0.5, 0)), plwd = 2, plty = 1,
    pfcol = scales::alpha(color, c(0, 0.5, 0, 0)), pcol= color, plty = c(2, 1, 0, 0), plwd = 2, pty = c(32, 16, 16, 16),
    # Customize the grid
    cglcol = "grey", cglty = 1, cglwd = 0.8,
    # Customize the axis
    axislabcol = "black", 
    # Variable labels
    vlcex = vlcex, vlabels = gsub(" ", "\n", vlabels),
    caxislabels = caxislabels, calcex = 0.6,
    title = title
  )
}

create_sequence <- function(min_val, max_val, n) {
  seq(min_val, max_val, length.out = n)
}

for (i in 1:dim(speciesPara)[1]) {
  png(paste0(rownames(speciesPara)[i], ".png"), width=1400, height=1200, res=300)
  #png(paste0("test", ".png"), width=1400, height=1200, res=300)
  max_min <- rbind(rep(max(speciesPara[i,], na.rm = T), dim(speciesPara)[2]), rep(min(speciesPara[i,], na.rm = T), dim(speciesPara)[2]))
  rownames(max_min) <- c("Max", "Min")
  colnames(max_min) <- colnames(speciesPara[i,])
  mean_value <- rowMeans(speciesPara[i,], na.rm=T)
  speciesPara_below <- speciesPara[i,]; speciesPara_below[speciesPara_below > mean_value] <- NA
  speciesPara_above <- speciesPara[i,]; speciesPara_above[speciesPara_above < mean_value] <- NA
  data_each_Para <- rbind(max_min, rep(mean_value, length(speciesPara[i,])), speciesPara[i,], speciesPara_below, speciesPara_above)
  # Reduce plot margin using par()
  op <- par(mar = c(2, 2, 2, 1))
  # title_name <- data.frame(x=1, y=1, label=rownames(speciesPara)[i])
  create_beautiful_radarchart(data_each_Para, color = c("grey","#00AFBB", "#6db3fc", "#f99690"),
                              caxislabels = create_sequence(min(speciesPara[i,], na.rm = T), max(speciesPara[i,], na.rm = T), 5),
                              title = NA)
  #title(rownames(speciesPara)[i], cex.main=0.8)
  title(paraDef[i], cex.main=0.8, line = 0.75) # adjust size and location
  vp <- viewport(width = 1, height = 1.15, x = 0.521, y = 0.5)
  donut_plot <- ggplot(donut_data, aes(ymax=ymax, ymin=ymin, xmax=3.75, xmin=3, fill=category)) +
    geom_rect(alpha = 0.4) +
    geom_text(aes(x = 3.85, y = label_position, label = category, angle = mid_angle), size = 2.5, color = "grey30") +  # Add category labels
    #geom_text(data=title_name, aes(x=1, y=1, label=label), inherit.aes = FALSE, size=2.5) +
    scale_fill_brewer(palette=8) +
    coord_polar(theta="y", start = 0.043*pi) +
    xlim(c(1, 4)) +
    theme_void() +
    theme(legend.position = "none")
  donut_plot
  
  # Plot ggplot2
  print(donut_plot, vp = vp, newpage = FALSE)
  
  par(op)
  dev.off()
}

png_names <- list.files(pattern = "*.png")  # Adjust accordingly
image_list <- lapply(png_names, image_read)
animation <- image_join(image_list)
animation <- image_animate(animation, fps=0.5)  # Adjust fps (frames per second) as needed
image_write(animation, "parameters_radarplot.gif")


## cluster analysis
set.seed(123)
data_para <- data.frame(t(speciesPara))
data_category <- list(data_para[,5:13], data_para[,c(5,6,8,10,11)], data_para[,c(9,13)], 
                      data_para[,c(5,6,8,9,10,11,13)], data_para[,c(1,3,4,5)], data_para[,c(2,7,12)])
data_category_name <- c("Species", "Heat", "Water", "Heat&Water", "Phenology", "Growth")
#data <- data[!(rownames(data) %in% c("Watermelon", "African nightshade")), ]
for (i in 1:length(data_category)) {
  data <- data_category[[i]]
  means <- apply(data,2,mean)
  sds <- apply(data,2,sd)
  nor <- scale(data,center=means,scale=sds) # Normalization
  noc <- 6 # set the number of clusters
  
  # Hierarchical agglomerative clustering 
  # png(paste0("Hierarchical_clustering_",data_category_name[i],".jpeg"), width=1400, height=1200, res=300)
  dist_matrix <- dist(nor) # Calculate distance matrix
  hc <- hclust(dist_matrix, method="complete")
  #hc$labels <- rep("", length(hc$labels))
  clusters <- cutree(hc, k=noc)
  plot(hc, labels = rownames(data), hang = -1, cex = 0.7,
       main = paste0("Cluster Dendrogram - ", data_category_name[i]), cex.main = 1,
       xlab = "", ylab = "Distance", ann = FALSE)
  title(paste0("Cluster Dendrogram - ", data_category_name[i]), cex.main = 1)
  rect.hclust(hc, k=noc)  # Draws rectangles around the clusters
  dev.off()
  
  # png("PCA.jpeg", width=1400, height=1200, res=300)
  # pca <- prcomp(nor)
  # plot(pca$x[,1:2], col=clusters, pch=16, 
  #      xlab = "Principal Component 1", ylab = "Principal Component 2",
  #      main="2D PCA colored by clusters")
  # legend("bottomleft", legend=paste("Cluster", 1:noc), fill=1:noc, cex = 0.7)
  # dev.off()
  
  # K-means clustering
  png(paste0("Kmeans_clustering_",data_category_name[i],".jpeg"), width=1800, height=1800, res=300)
  kmeans_result <- kmeans(nor, centers=3)
  datadistshortset<-dist(nor,method = "euclidean")
  hc1 <- hclust(datadistshortset, method = "complete" )
  pamvshortset <- pam(datadistshortset, noc, diss = FALSE)
  clusplot(pamvshortset, shade = FALSE,labels=2,col.clus="blue",col.p="red",span=FALSE,main=paste0("Cluster Mapping - ", data_category_name[i]),cex=0.8)
  dev.off()
}

dist_matrix <- as.matrix(dist(nor), diag = FALSE, upper = FALSE)

# Set the upper triangle and diagonal to NA
# dist_matrix[!lower.tri(dist_matrix, diag = FALSE)] <- NA  # Set diag = TRUE if you want the diagonal as NA too
diag(dist_matrix) <- NA

# Melt the matrix into a long data format
matrix_long <- melt(dist_matrix)
matrix_long$ValueFactor <- factor(cut(matrix_long$value, breaks=7, labels=FALSE), exclude = NA)  # Convert to discrete factor
matrix_long$ValueFactor <- droplevels(matrix_long$ValueFactor)
matrix_long$ValueFactor <- factor(matrix_long$ValueFactor, exclude = "NA")

# Define a function to create a jet-like color palette
jet.colors <- function(n) {
  colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))(n)
  colors <- colorRampPalette(c("#318AD0","#345F5D","#47A664","#FEDE24","#E66A1F","#F74EA7","#8872B3"))(n)
  colors <- colorRampPalette(c("#149580", "#23BB98", "#52D5B2", "#82DFC6", "#B3E7D0", "#FFEE00", "#FBB806", "#F6830C", "#F24D11", "#ED1717"))(n)
  return(colors)
}

# Create the heatmap
p <- ggplot(matrix_long, aes(x=Var1, y=Var2, fill=ValueFactor)) +
  geom_tile() +  # use geom_tile() to create the heatmap tiles
  # scale_fill_distiller(palette = "Spectral", direction = -1, name = "Distance", na.value = "transparent") +  # Reverse the Spectral palette
  scale_fill_manual(values = c("1" = "#149580","2" = "#23BB98","3" = "#52D5B2","4" = "#FFEE00", "5" = "#FBB806", "6" = "#F6830C","7" = "#ED1717"), 
                    labels = c("1", "2", "3", "4", "5", "6", "7", ""), 
                    name = "Distance", na.value = "transparent") +
  #scale_fill_gradientn(colours = brewer.pal(11, "Spectral", direction = -1), name = "Distance", na.value = "transparent") +
  labs(x = "", y = "", fill = "Value") +
  theme_classic() +
  ggtitle("Crop Distance") + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))  # Rotate x-axis labels

ggsave(filename = paste0("/Users/meijian/Work/SIMPLE20181102/parameters_radarplot/distanceHeatmap.png"), plot = p, width = 7, height = 6, units = "in", dpi = 300)

